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AN  ITERATIVE  SOLUTION  TO  THE  GENERALIZED  EIGENVALUE  -  EIGENVECTOR  PROBLEM 

I.   Statement  of  the  Problem 

The  "generalized  eigenvalue  problem"  may  be  stated  as  follows: 
Find  those  values  of  the  complex  parameter,  z,  which  will  provide 

a  non-trivial  solution  to  the  system  of  linear  equations, 


A11(z)X1  +  A12(z)X2  +  ...  +  Aln(z)Xn  =  0 


A21(z)X1  +  A22(z)X2  +  ...  +  A2n(z)Xn  =  0 


(1) 


A  , (z)X.  +  A  _(z)X9  +  ...  +  A  (z)Xn   =  0, 
nl    1    nz    I  nn    n 


This  can  be  written  more  concisely  as 

(2)  A(z)X  =  0 

where  A  is  an  n  x  n  matrix,  the  elements  of  which  are  arbitrary  complex 
functions  of  z,  and  X  is  a  complex  column  vector  of  n  elements.   The  set 
£z.:A(z.)X.  =  Oj  is  the  set  of  eigenvalues  and  the  X.  are  called  the 
associated  eigenvectors. 

The  special  case  where  the  A..(z)'s  are  polynomials  in  z  of  degree 
2  or  less  is  encountered  most  often  in  practical  work.   The  rest  of  the 
discussion  in  this  paper  will  be  confined  to  problems  of  this  type.   However, 
it  should  be  pointed  out  that  the  general  method  described  here  is  not  so 
restricted. 


Under  the  above  restriction  we  can  rewrite  A(z)  as 

(3)  A(z)  =  Aqz2  +  A  z  +  A2 

where  the  A.'s  are  constant  complex  matrices.   We  wilj  need  this  form 
in  the  later  discussion  on  scaling. 

II.   Method  of  Solution 

The  heart  of  the  method  described  here  is  the  iterative  root-finding 
technique.   After  a  literature  search,  in  which  much  background  information 
was  gained  from  Wilkinson  (1),  and  some  experimentation,  it  was  decided  that 
a  method  developed  by  J.  F.  Traub  (2)  would  be  used  in  this  code.   Muller's 
method  (7)  has  been  mentioned  by  a  number  of  authors  in  connection  with 
iterative  solutions  in  the  complex  plane;  however,  Traub' s  technique  was 
chosen  because  of  its  simplicity  and  ease  of  computation.   The  methods  are 
similar;  both  are  of  second  order  and  have  the  same  asymptotic  error  constant. 
The  main  difference  is  that  Traub  uses  the  Newtonian  form  of  the  interpolating 
polynomial  and  Muller  uses  the  Lagrangian  form. 

In  essence  the  technique  is  to  perform  a  quadratic  fit  using  the  last 
three  iterants  and  their  functional  values.  This  algorithm  is  defined  as 
follows: 

(4)  z...  =  z.  -2f./(w+  (w2-4f.  f[  z.  z.  .  z.  .]}  1/2) 

l+l    l    l   —  *»     l.i,  l-l,  i-2  "* 

where 

(5)  w=f[  z.  z.   ]   +  (z.-z.   )f[  z.   z.    z.   ] 

l,  l-l  .      l   l-l     l,   i-l,  1-2 

and 

(6)  f[z   z.  .  z.  .]  s  (f[z.  z.  .]  -  f[z.  .  z.  J)/(z.-z.  _) 

i,  i-l,  i-2        l,  i-l       i-1j  i-2     l   i-2 

The  variable  sign  in  equation  (4)  is  chosen  so  that  the  denominator  has 
the  larger  magnitude.  _, 
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Convergence   is   determined  by  either  of   the   two   tests, 
(7)  lzi-Zi-l|/|zil    <el 


or 

(8)  lfil<e2      ' 

If  either  of  these  inequalities  is  met  the  process  is  considered  to 
have  converged  to  a  root. 

In  order  to  avoid  convergence  to  a  previously-calculated  root,  a 
method  of  suppressing  zeros,  first  suggested  by  G.  E.  Forsythe  (Reference  3, 
page  164)  was  incorporated  into  the  code.   The  first  root,  z-      is  found 
as  a  zero  of  f(z)  and  the  remaining  roots,  z.,  are  found  as  zeros  of 
g. (z) ,  where 

(9)  g.(z)  s  f(z)/TT  (z-z.)  ,    i>  1   . 

1  k=l     * 

One  of  the  problems  encountered  in  a  program  of  this  type  is  floating- 
point exponential  overflow  when  working  with  matrices  of  large  order.   To 
alleviate  this  problem  an  automatic  scaling  feature  was  devised.   This 
scaling  is  done  in  two  parts,  a_   priori  and  a_   posteriori.   The  a^  priori 
scaling  is  done  at  the  beginning  of  the  problem  and  is  carried  out  in  the 
following  manner.   The  coefficients  of  each  polynomial  element  in  a  row  of 
the  matrix  A  are  summed  and-a  maximum  is  found.   In  a  binary  computer  this 
maximum  may  be  expressed  exactly  as  bx2  (.   Then  each  coefficient  in  every 
polynomial  in  that  row  is  divided  by  2  .   This  operation  is  carried  out 
for  each  row  in  the  matrix. 

A  posteriori  scaling  is  done  only  when  an  overflow  occurs.   If,  at 
the  time  of  the  overflow,  the  estimate  to  the  root  being  calculated  is  some 
number,  say  cx2  ,  then  A(z)  (defined  in  equation  (3)   )  is  replaced  by 


(10) 


A'(z')  =  A  z  +  —  z  +  — — _ 
2k     (2k)2 


and  a  root,  z',  is  then  found  for  the  new  problem.   The  root  of  the 

original  problem  is  then  calculated  as 

(11)  z.  =  2kz* 
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The  scaled  problem  is  started  by  taking  the  last  three  iterants  in  the 
original  problem  and  dividing  them  by  2  .   After  a  root  is  found  the 
matrix  is  returned  to  its  original  form  and  the  program  proceeds  to  find 
the  next  root. 

In  many  cases  it  is  desirable  to  obtain  the  eigenvectors  associated  with 
the  eigenvalues.   An  option  to  perform  these  calculations  is  available 
in  the  program.   Because  there  exists  no  procedure  which  will  absolutely 
assure  the  finding  of  an  eigenvector,  the  code  will  make  2n  attempts, 
where  n  is  the  order  of  the  matrix,  at  computing  an  eigenvector  for  each 
eigenvalue.   This  procedure  is  carried  out  as  outlined  below. 

An  assumption  is  made  that  the  elements  of  the  vector  may  be  written 
as  multiples  of  a  single  element  in  the  vector.   One  possible  representation 
of  this  is 


(12) 


X.  =  (x    b  x   b_x    ...,  b   X  .) 

l     li,   1  2i,  2  3i,      n-1  ni 


Using  this  assumption  the  code  makes  its  first  attempt  by  setting  X..  .  =  1 
and  solving  the  system  of  equations  (13): 

/~a2l\ 


\ 


(13) 


a32  a33*  '  *  a3n 


a   a   ...  a 
n2  n3       nn 


X 


3i 


/  \w 


where  a.,  is  the  function  A.,  evaluated  at  z.. 

The  expression,  A(z)  X,  is  then  calculated  using  this  vector 
and  the  result  compared  to  the  zero  vector.   If  the  vector,  X,  is 
found  to  be  unacceptable  a  second  attempt  is  made  by  replacing  a„  by  a... 
k=l,2,...,n,  in  (13).   If  this  fails  then  x   is  set  equal  to  unity 


and  the  system  (14)  is  solved, 


/ 


all  a13' 


(14) 


-\  I %"\ 


/-< 


a31  a33  *  '  *a3n 


3i 


a   a  v 
nl  n3 


.  a 


nn 


/  UJ 


121 


-a 


32 
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and  checked  as  indicated  above.   A  fourth  attempt  may  be  made  by  replacing, 
in  (14),  a   by  a   for  k=l,2,...,n.   The  process  continues  for  all  X..  or 

•J  K.  Z  rC  J  1 

until  an  eigenvector  is  found.   If  no  eigenvector  is  found  a  statement  to 
that  effect  is  printed  along  with  the  results  of  the  final  attempt. 


III.   General  Discussion  of  Method 

A  code,  using  the  method  described  in  this  paper,  has  been  written 
in  FORTRAN-63  for  use  on  a  CDC  1604.   This  program  has  been  checked  out 
on  a  variety  of  problems,  several  of  which  are  shown  in  the  Section  IV. 
The  accuracy  of  the  process  has  been  shown  to  be  quite  good.   This  may 
be  attributed  to  several  factors.   First,  only  the  original  input  matrix 
is  used  for  the  evaluation  of  the  function  each  time,  and,  second,  all 
scaling  is  done  by  powers  of  two  (  Reference  (8) ,  page  98  ) .   The  method 
has  been  shown  to  be  reasonably  fast,  although  a  good  relative  comparison 
is  impossible  because  the  number  of  iterations  required  is  a  function  of 


the  estimate  provided.   Unfortunately  this  information  generally  has  not 
been  included  in  previous  literature  on  this  subject.   Because  the  number 
of  computations  necessary  to  determine  the  next  estimate  to  the  root  is 
so  small  compared  to  those  required  to   evaluate  the  determinant,  the  time 
required  for  each  iteration  can  be  considered  substantially  equal  to  the 
time  necessary  for  the  latter.   The  speed  of  the  program  is  further  enhanced 
by  the  fact  that  through  the  use  of  good  estimates  only  eigenvalues  of 
interest  need  be  computed,  although  convergence  seems  to  result  from  almost 
arbitrary  estimates.   Also  because  of  the  feature  which  forces  avoidance  of 
previously  found  zeros,  it  is  possible  to  restart  the  program  and  by-pass 
any  such  prior  calculations. 

IV.   Numerical  Examples 

Most  of  the  examples  used  for  this  section  have  been  taken  from  papers 
on  the  standard  eigenvalue  problem.   This  was  necessary  because  of  the 
lack  of  comparative  data  on  the  generalized  problem.   However,  these  examples 
should  be  adequate  to  demonstrate  the  ability  of  this  method  to  find  path- 
ological roots  and  also  to  get  some  comparison  as  to  its  accuracy.   All 
problems,  using  the  method  described  in  this  report,  were  run  on  a  CDC  1604 
which  has  a  48-bit  word  length. 


Example  1 

(See  Parlett  (5)   ).   These  results  were  obtained  using  Laguerre's 
method  on  an  IBM  7090  which  has  a  36-bit  word.   A  is  defined  as 
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-49 

29  \ 

43 
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52 

44 

59 

-23 
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208 

208 

411-z 

208 

208 

208 

99-z 

-911 

208 

-911 

99-z/ 

Table  of  Eigenvalues 


v.;--- 

Actual  (8  digits) 

Parlett  (5) 

Present  Method 

1020.0490 

1020.0500 

1020.0490 

1020.0000 

1019.9997 

1020.0000 

1019.9020 

1019.9019 

1019.9019 

1000.0000 

1000.0001 

1000.0000 

1000.0000 

999.99999 

1000.0000 

.098048640 

.098045509 

0.098048646 

0.0 

-0.0000094 

-0.000000003 

-1020.0490 

-1020.0490 

-1020.0490 

Total  running  time:   16  seconds 

Average  number  of  iterations  per  eigenvalue: 


10.5 


Example  2 


(See  Barlow  and  Jones  (4)   ) .   Their  work  was  done  on  an  IBM  7040 
(36-bit  word)  using  an  extension  of  the  classical  secant  method  to  the 
complex  field.   This  paper  (4)  also  gives  results  obtained  by  Eberlein  (6) 
using  a  Jacobi-like  method. 

6 

9 

6-z 

5 

5 


A  = 


15-z 

11 

1 

3-z 

7 

6 

7 

7 

17 

12 

-9 

-15 

—3 

-8 

-3 

-11 

-3-z 

-11 

-10 

-16-z 

Table  of  Eigenvalues 


Actual  (6) 

Eberlein  (6)  1 

Barlow  &  Jones  (4) 

Present  Method 

-i-:o 

-.999 

-1.0000004 

-.999999999 

1.50016  +  3.57064i 

1.505  +  3.571 

1.49850  +  3.57027i 

1.50001  +  3.57073i 

1.50016  +  3.57064i 

1.495  +  3.57i 

1.50150  +  3.571101 

1.50003  +  3.570691 

1.50016  -  3.570641 

1.505  -  3.57i 

1.49996  -  3.571271 

1.50001  -  3.57074i 

1.50016  -  3.570641 

1.495  -  3.57i 

1.49864  -  3.57023i 

1.49999  -  3.57070i 

Total  running  time:   4.38  seconds 

Average  number  of  iterations  per  eigenvalue:   6 


Example  3 

(See  Parlett  (5)   ).   This  example  was  used  to  demonstrate  the 
ability  of  the  code  to  extract  only  a  few  desired  roots. 


/  C    2C\ 

\4C    3CJ 


A  =  (B  -  zl) 

5D  -  d' 


C  = 


(S"S) 


D  = 


The  true  eigenvalues  of  this  system  (Eberlein  (6)   )  are  60+20i, 
45+15i,  30+10i,  15+51,  -12+41,  -9+3i,  -6+2i,  -3+i. 


Initial  Guess 

Eigenvalue  Found 

-2.5  -  0.6251 
-2.5  +  0.6251 
-5.0  -  1.25i 
-5.0  +  1.251 

-3.0000006  -  9.9999947i 
-3.0000005  +  l.OOOOOOOi 
-6.0000055  -  1.9999948i 
-5.9999975  +  2.0000003i 

Total  running  time:   55.28  seconds 

Average  number  of  iterations  per  eigenvalue:   17.5 


Example  4 


In  this  example  both  the  eigenvalues  and  eigenvectors  were  found, 


A  = 


1-z       -2  3  -2 

1  5-z  -1  -1 

2  3  2-z  -2 
2-2  6  -3-z 


Table  of  Eigenvalues  and  Eigenvectors 


True  Values 

Computed  Values 

Eigenvalues 

Eigenvectors 

Eigenvalues 

Eigenvectors 

-1 

1 
0 
0 

1 

-1.00000000 

1.0 

0.0 
0.0 
1.0 

2 

1 
-6 
-9 
-8 

1.999939 

1.0 
-5.9993 
-8.9991 
-7.9993 

2 

1 
-6 
-9 
-8 

2.00007 

1.0 
-6.0007 
-9.001 
-8.0008 

2 

1 
-6 
-9 
-8 

1.99998 

1.0 
-5.9998 
-8.9998 
-7.9998 

Total  running  time:   3.83  seconds 

Average  number  of  iterations  per  eigenvalue:   8.25 
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Example  5 


This  is  an  example  of  the  generalized  problem. 


2z  +14z-88    -z  -13z-22 


z  +2z-99 


A  = 


-3z  +21z-36     6z  -6z-36 


4z  -48z+108 


8z  +8z-160     z  +7z+10     -z  +4z+45 


Table  of  Eigenvalues  and  Eigenvectors 


True 

Values 

Computed  Values 

Eigenvalues 

Eigenvectors 

Eigenvalues 

Eigenvectors 

-11 

1 
-29/6 
153/40 

-10.99999999 

: 

1.0 
-4.833333333 
3.82500000 

-5 

10/33 
1 
-27/154 

-5.0 

0.3030303030 
1.0 
-1.753246753 

-2 

0 
1 
0 

-2.0 

0.0 
1.0 
0.0 

3 

0 

1 

-5/6 

3.0 

-0.0000000001 

1.0 
-0.8333333333 

4 

1 
0 
0 

4.0 

1.0 
0.0 
0.0 

9 

0 
0 

1 

9.0 

0.0 
0.0 
1.0 

Total  running  time:   3.85  seconds 

Average  number  of  iterations  per  eigenvalue: 


4.0 
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Example  6 

This  final  example  of  the  generalized  problem  is  constructed  so 
that  it  can  also  be  solved  as  a  standard  problem.   It  is  defined  as  follows , 

i.)    As  (B+zC)  =  0, 
where  B  and  C  are  12x12  matrices  with  constant  elements.   By  multiplying 
on  the  left  by  C   a  transformation  is  made  to  the  standard  problem, 

ii.)   A'  s  (BC^+zI)  =  0. 
This  problem  has  been  solved  in  both  forms  by  the  technique  discussed  here. 
In  addition,  it  was  solved  using  the  algorithm  in  (5),  so  that  the  results 
might  be  compared. 

The  matrices  B  and  C  are  defined  on  the  following  page. 


3*$ 

12 


The  matrices  B  and  C  were  defined  as: 


(  10 

-15 

1 

-2 

-81 

91 

-69 

14 

62 

-36 

-20 

99> 

22 

-46 

-25 

-85 

-3 

89 

27 

53 

-93 

-34 

52 

19 

24 

48 

-22 

-97 

-76 

-64 

-15 

24 

49 

-32 

30 

-19 

-42 

-93 

-6 

61 

7 

16 

3£ 

-53 

-71 

57 

0 

-74 

37 

-39 

-81 

16 

-6 

91 

6(1 

-81 

49 

60 

14 

-6 

-77 

-6 

11 

-42 

27 

53 

18 

-70 

-90 

-15 

21 

-81 

99 

-72 

56 

69 

98 

31 

71 

18 

-44 

48 

-63 

-21 

-96 

-91 

-5 

7 

18 

20 

-94 

-56 

69 

-60 

-18 

-84 

-89 

14 

-63 

-10 

-17 

-18 

57 

84 

-25 

12 

58 

-44 

-85 

-36 

53 

53 

53 

-59 

-38 

62 

8 

-17 

-16 

11 

28 

69 

-88 

33 

-70 

79 

-56 

-5 

90 

-31 

-1 

85 

^-63 

-40 

-48 

-3 

49 

-69 

-18 

-72 

52 

-20 

12 

-90, 

'-81 

-:72 

-4 

-96 

24 

-82 

66 

14 

-76 

14 

13 

8l\ 

29 

-20 

68 

26 

-46 

-20 

89 

81 

-86 

-12 

-92 

57 

0 

57 

39 

66 

-84 

-40 

32 

61 

-98 

-96 

64 

64 

5 

-4 

-25 

26 

-44 

44 

-37 

63 

45 

66 

75 

66 

-91 

26 

-64 

-94 

26 

25 

39 

-22 

71 

64 

91 

42 

0 

-4 

-87 

-77 

42 

35 

-74 

-99 

-81 

-42 

43 

-76 

0 

69 

-62 

56 

86 

88 

76 

36 

-84 

-93 

76 

-65 
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These  three  cases  were  run  with  a  relative  error  bound  of  10  "  and 
the  results  from  all  of  them  were  identical  to  six  or  seven  significant 
figures.   The  eigenvalues  found  are  listed  in  the  table  below. 


Table  of  Eigenvalues 


.956212 
■2.67567 
■  .0213507 

•  .466901 

•  .626829  +  2.56009 

•  .274574  +  .327271 
.254714  +  .457145 

•  .796302  +  1.04215 


Total  time  for  present  method:   104  seconds 
Average  number  of  iterations  per  eigenvalue:   17.2 
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V.   Areas  of  Further  Investigation 

Some  additional  features  are  being  considered  for  inclusion  in  the 
present  code.   However,  at  the  time  of  this  report  all  problems  involving 
these  revisions  had  not  been  resolved. 

It  is  felt  that  some  process  which  could  quickly  determine  the 
approximate  location  of  eigenvalues  may  be  useful,  especially  when  large 
order  problems  are  undertaken.   It  has  not,  as  yet,  been  determined  what 
method  might  be  most  profitable  for  this  purpose. 

I      An  automatic  entry  into  double  precision  would  be  desirable  for  parts 
of  the  computations.   But  it  will  be  necessary  to  devise  some  general  test 
to  indicate  when  this  move  could  best  be  made  without  wasting  computer 
time. 

»A  test  to  find  when  the  domain  of  indeterminancy  has  been  reached  in 
the  search  for  a  root  would  be  a  valuable  addition   to  those  existing 
tests  for  acceptance  of  a  root. 
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Calculate  and 
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APPENDIX  III 


IDENTIFICATION; 

Title:       A  Method  for  the  Solution  of  the  "Generalized 

Eigenvalue  Problem"  with* Polynomial  Matrix  Elements 

CO-OP  ID:  F4-NPGS-GENEIG        (1604  F-63) 

Category:  Mathematical  Subroutine 

Programmer:  R.  D.  Brunell 

Date:  June  10,  1966 

PURPOSE: 

To  find  those  values  of  the  complex  parameter,  z,  which 
satisfy  the  system  of  equations: 

A11(z)x.  +  A1.(z)x.  +  ...  +  A_  (z)x  =  0 
11     1     12    2  In    n 

A21(Z)X1  +  A22(Z)X2  +  '••  +  A2n(z)Xn  =  ° 


(1) 

Anl(2)xi+An2(z)x2+  -  +  Ann(z)Xn  =  °- 


where  the  a..'s  are  complex  polynomials  of  degree  2  in  z, 
This  may  be  written  more  concisely  as 


H(z)X  =  0, 

where  H(z)  is  an  n  x  n  matrix  whose  elements  are  the  polynomials 
in  z,  and  X  is  a  complex  column  vector  of  n  elements.   The  elements 
of  the  set  £z.:H(z.)X.  =  0J  are  called  the  eigenvalues,  and  the 
X.  are  called  the  associated  eigenvectors. 

The  matrix  H(z),  because  of  the  form  of  its  elements,  may  be 
wsitten  in  the  following  form: 

H(z)  =  A  z2  +  A  z  +  A    , 

where  the  n  x  n  matrices  A.  have  constant  complex  elements. 

Equation  (3)  will  be  useful  later  when  discussing  the  input  for  the 
subroutine. 
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C.   USAGE: 


(4) 


1.  Set  up  for  GENEIG 

GENEIG  is  a  subroutine  and  the  user  must  supply  a  main 

program  which  will  set  up  the  input  data  for  GENEIG.   All 

communication  to  and  from  this  subroutine  is  done  through 

common  storage,  and  the  following  DIMENSION,  COMMON  and  TYPE 
statements  must  be  declared  by  the  user. 

DIMENSION    A(20,20,3) ,EIGG(40,2) ,EVAL(40) ,EVEC(20,40) 

COMMON       A,EIGG,EVAL,EVEC,NN,IORD,IE, 

NEIG , IMAX, ISAC , IPRINT , IVEC , EPI , EP2 , EP3 

TYPE  COMPLEX  A,EIGG,EVAL,EVEC 

2.  Input  for  GENEIG 

The  three-dimensional  array,  A,  must  contain  the  elements 
of  the  matrix  H  in  the  manner  described  below: 


a0).  .    =  A(I,J,1),  iH  .  =  A(I,J,2)  and  a2\  .    =  A(I,J,3), 
ij  ij  ij 

(k) 
where  a  . .  is  the  ii-th  element  of  A,  . 
ij  k 

The  two-dimensional  array,  EIGG,  must  contain  the  estimated 
value  for  the  eigenvalues  as  follows: 

EIGG(I,1)  will  be  considered  as  the  estimate  to  the  value  of  the 
i-th  eigenvalue,  and  EIGG(I,2)  will  be  considered  as  the  probable 
precentage  error  in  that  guess.   The  first  three  values  necessary 
to  start  the  process  for  the  i-th  value  are  computed  by 

Jr.    =   EIGG(I,1)  -  EIGG(I,2)*EIGG(I,1) 

(2) 
z  .  =  EIGG (1,1) 

(3) 

z  .  =  EIGG(I,1)  +  EIGG(I,2)*EIGG(I,1) 

The  array  EVEC  will  contain  the  associated  eigenvectors.   The 
eigenvector  associated  with  the  i-th  eigenvalue  will  be  found  in 
EVEC(J,I),  where  J  =  1,2,...,NN. 
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NN  is  the  order  of  the  coefficient  matrices. 

IORD  is  the  degree  of  the  polynomial  elements. 

IE  is  the  number  of  the  first  eigenvalue  to  be  found. 

NEIG  is  the  number  of  the  last  eigenvalue  to  be  found. 

IMAX  is  the  maximum  number  of  iterations  to  be  allowed. 

ISAC  is  a  flag  for  initial  scaling.   If  ISAC  =  0  no  scaling 
will  occur. 

IPRINT  is  a  flag  for  printing  at  every  iteration.   If  IPRINT  =  0 
only  the  final  results  will  be  printed. 

IVEC  is  a  flag  for  the  eigenvector  option.   If  IVEC  =  0 
no  eigenvectors  will  be  computed. 

EPI  is  the  relative  error  bound  on  the  eigenvalues.   If 

|z .  -  z._1|  /  |z .|  ^ EPI  the  problem  is  considered  to  have  converged, 

-4  -8 

A  value  of  10   ^  EPI  £  10   has  been  found  to  be  satisfactory 

for  most  cases. 

EP2  is  an  absolute  error  bound  on  the  function.   If  JH(z.)/<dEP2 

the  problem  is  considered  to  have  converged.   A  value  of  EP2  fti 10 
appears  to  be  satisfactory  for  those  cases  for  which  this  test  becomes 
effective. 

EP3  is  a  relative  error  bound  on  the  eigenvector.   If  B  is  the  maximum 

element  of  H(z.),  then  if  JH(z.)X.|  <  |eP3*b|  the  eigenvector 

-4 
is  considered  acceptable.   A  value  of  EP3  £i  10   is  satisfactory. 

3.   Calling  Sequence: 

Once  the  common  storage  has  been  set  up  in  the  user's  program 
it  is  only  necessary  to  insert  the  instruction  CALL  GENEIG.   The 
subroutine  will  then  compute  all  eigenvalues  and  eigenvectors,  that 
have  been  indicated  by  the  input,  before  returning  control  to  the 
calling  program. 


■15 


1 

i 
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■ 


Space  Required:   4217..  n  excluding  common  storage 


Common  Storage  Required:   4251 


10 


Print-Outs:   If  IPRINT  =  0  then  at  each  iteration  the  number  of 

the  iteration,  the  present  value  of  z.  and   H(z.)   will  be  printed. 

When  an  eigenvalue  is  accepted,  a  message  to  this  effect  will  be 
printed  along  with  the  accepted  value  and  the  function  value. 

If  IVEC  =  0  and  an  eigenvector  is  found  it  will   be  printed. 

However,  if  no  acceptable  eigenvector  is  found,  a  statement  to 

this  effect  will  be  printed  along  with  the  value  of  the  final  attempt 

Error  Stops:   If  IORD  is  found  to  be  other  than  an  integer  1  or  2 
the  program  will  print  the  value  of  IORD  and  stop  (Exit  to  Monitor) . 
If  the  determinant  evaluator  is  unable  to  complete  its  operation  the 
code  will  indicate  that  this  has  happened  and  stop  (Exit  to  Monitor) . 

Input  and  Output  Formats:   Not  applicable. 

Input  and  Output  Tape  Mountings:   Output  is  to  Logical  unit  51 
(Standard  on  CO-OP  Monitor). 


10.  Selective  Jump  and  Stop  Settings:   Not  applicable. 

11.  Timing  Examples: 


Case 

1 

2 

3 

,   4 

5 

,  6 

Size  of  Matrix 
No.  of  Eigenvalues 
No.  of  Eigenvectors 
Average  No.  of  Iterations 
Total  Time  (seconds) 

3x3 

6 

6 
4.0 
3.85 

4x4 

4 

4 
8.25 
3.85 

5x5 

5 

0 
6.0 
4.38 

8x8 

8 

0 

10.5 

16.0 

» 

12x12 
12 
0 
18.7 
104. 

16x16 

4 

0 

17.5 

*  55.28 

1 

12.   Accuracy: 

Results  from  this  code  have  been  compared  with  results  taken  from 
several  papers  on  the  standard  eigenvalue  problem  and  they  have  been 
as  accurate  or  more  accurate  on  all  tests.   Also,  several  problems 
were  run  for  which  exact  results  were  known.   On  problems  with  well- 
defined  roots  nine  or  ten  figure  accuracy  was  observed.   When  pathological 
roots  were  encountered  the  accuracy  of  this  code  compared  very  favorably 
with  other  methods  tested. 
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13.  Cautions  to  User:   The  DIMENSION,  COMMON  and  TYPE  statements 
as  shown  in  Section  C  must  be  used  exactly  as  indicated. 

14.  Equipment  Configuration:   Standard  CO-OP  Computer  with  FORTRAN  63 
Compiler  and  CO-OP  Monitor  System. 

15.  References: 

(a)  R.  D.  Brunei 1,  An  Iterative  Solution  to  the  Generalized 

Eigenvalue-Eigenvector  Problem,  Technical  Report/Research 
Paper  No.  69,  U.S.  Naval  Postgraduate  School,  Monterey, 
California,  1966. 

(b)  J.  F.  Traub,  Iterative  Methods  for  the  Solution  of  Equations, 

Prentice-Hall,  Inc.,  Englewood  Cliffs,  New  Jersey,  1964. 

(c)  I.  Tarnove,  "Determination  of  Eigenvalues  of  Matrices  Having 

Polynomial  Elements,"  Journal  of  SIAM,  June  1958. 


D.   METHOD: 

Only  a  general  discussion  of  the  method  will  be  presented  here. 
For  a  more  detailed  description  see  Reference  (a). 

The  heart  of  the  method  used  in  this  subroutine  is  the  iterative 
root-finding  technique,  which  was  developed  by  J.  F.  Traub  (b) .   In 
essence  the  technique  is  to  perform  a  quadratic  fit  using  the  last 
three  iterants  and  their  functional  values.   The  new  estimate  to  the 
root  is  defined  as  that  root  of  the  fitting  quadratic  which  is  closer 
to  the  last  iterant.   The  process  converges  with  almost  arbitrary 
estimates;  however,  through  the  use  of  good  estimates  a  considerable 
saving  in  computer  time  can  be  realized.   Fewer  iterations  will  be 
necessary  and  only  those  eigenvalues  of  interest  need  be  computed. 

In  order  to  avoid  convergence  to  a  previously-calculated  root 
a  method  of  suppressing  zeros,  first  suggested  by  G.  E.  Forsythe 
(Reference  (c)  page  164) ,  was  incorporated  into  the  code.   This 
feature  can  also  be  used  to  eliminate  re-calculation  of  roots  found 
on  a  previous  run,  especially  when  they  may  lie  pathologically  close 
to  a  desired  root.   This  can  be  accomplished  by  loading  the  m  known 
roots  into  the  first  m  positions  in  the  EVAL  array  and  setting  IE 
equal  to  m  +  1. 

To  alleviate  the  problem  of  floating-point  exponential  overflow 
a  two-part,  a  priori  and  a  posteriori,  automatic  scaling  procedure 
was  devised  for  this  code. 
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The  a_  priori  scaling  is  optional  and  is  controlled  by  the 
flag  ISAC.   This  scaling  occurs  only  once  and  is  done  prior  to 

any  calculations  for  eigenvalues.   The  coefficients  of  each  ! 

polynomial  element  are  summed  and  a  maximum  ,B,  is  found  in  each 
row.   The  scaling  is  carried  out  by  dividing  each  coefficient  of 

every  polynomial  in  a  row  by  2    ,  where  2  ^- B  ^ 2 


; 


The  a.  posteriori  scaling  is  not  optional   and  will  occur  only 
when  an  overflow  is  sensed.   At  this  point  the  problem  is  trans- 
formed so  that  the  scaled  eigenvalue  will  be  exactly  the  original 

g 

k+l         k         k+1  4 

eigenvalue  divided  by  2    ,  where  2  -^  z.  ^  2    and  z.  is  the  value  J 

i  i 

of  the  iterant  at  the  time  of  the  overflow. 

The  computation  of  the  associated  eigenvectors  is  an  optional 
feature  in  this  code.   The  basic  assumption  in  this  section  is  that 
the  eigenvalue  reduces  the  rank  of  the  original  matrix  by  1.   Under 
this  assumption  2NN  attempts  will  be  made  to  find  an  acceptable 
eigenvector.   If  all  attempts  fail,  a  printout  to  this  effect  will 
be  made. 


I 

i 

i 
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